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ABSTRACT 

We  consider  the  problem  of  feasibility  determination  in 
a  stochastic  setting.  In  particular,  we  wish  to  determine 
whether  a  system  belongs  to  a  given  set  F  based  on  a 
performance  measure  estimated  through  Monte  Carlo  sim¬ 
ulation.  Our  contribution  is  two-fold;  (i)  we  characterize 
fractional  allocations  that  are  asymptotically  optimal;  and 
(ii)  we  provide  an  easily  implementable  algorithm,  rooted 
in  stochastic  approximation  theory,  that  results  in  sampling 
allocations  that  provably  achieve  in  the  limit  the  same  perfor¬ 
mance  as  the  optimal  allocations.  The  finite-time  behavior 
of  the  algorithm  is  also  illustrated  on  two  small  examples. 

I  BACKGROUND 

We  consider  the  problem  of  feasibility  determination  in 
a  stochastic  setting.  In  particular,  we  wish  to  determine 
whether  a  system  belongs  to  a  given  set  F  based  on  a 
performance  measure  estimated  through  Monte  Carlo  sim¬ 
ulation.  While  it  is  an  interesting  problem  in  its  own  right, 
feasibility  determination  has  recently  attracted  much  atten¬ 
tion  from  researchers  in  ranking  and  selection  (R&S)  whose 
objective  is  to  select  the  best  system  in  the  presence  of  a 
stochastic  constraint. 

Ranking  and  selection  (R&S)  techniques  are  statistical 
methods  developed  to  select  the  best  system,  or  a  subset 
of  systems  from  among  a  set  of  alternative  system  de¬ 
signs.  R&S  via  simulation  is  particularly  appealing  as  it 
combines  modeling  flexibility  of  simulation  with  the  effi¬ 
ciency  of  statistical  techniques  for  effective  decision  making. 
Furthermore,  it  is  relatively  straightforward  to  satisfy  the 
underlying  technical  assumptions  of  these  techniques  in 
simulation  experiments,  which  also  allow  for  multi-stage 
sampling  as  required  by  some  R&S  methods. 

Due  to  randomness  in  output  data,  comparing  a  number 
of  simulated  systems  requires  care.  Suppose  we  conduct  n 
simulation  replications  for  each  of  r  designs.  Therefore,  we 
need  rn  simulation  replications.  Simulation  results  become 
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more  precise  as  n  increases.  If  the  precision  requirement 
is  high  in  is  not  small),  and  if  the  total  number  of  designs 
in  a  decision  problem  is  large  (r  is  large),  then  rn  can  be 
very  large,  which  may  easily  render  the  total  simulation 
cost  prohibitively  high  and  preclude  the  feasibility  of  using 
simulation  for  R&S  problems.  The  effective  employment 
of  the  simulation  budget  in  the  course  of  obtaining  a  good 
decision  is  therefore  crucial. 

The  rich  literature  on  R&S  has  been  developing  along 
two  principal  axes.  The  optimal  computing  budget  alloca¬ 
tion  (OCBA)  approach  has  its  roots  in  Dudewicz  and  Dalai 
(1975)  who  developed  a  two-stage  procedure  for  selecting 
the  best  design  or  a  design  that  is  very  close  to  the  best  sys¬ 
tem.  Rinott  (1978)  presents  an  alternative  way  to  compute 
the  required  number  of  simulation  replications  in  the  sec¬ 
ond  stage.  This  idea  has  been  ultimately  extended  to  fully 
sequential  algorithms.  Paulson  (1964),  Gupta  and  Pancha- 
pakesan  (1979),  Nelson  and  Matejcik  (1995),  Matejcik  and 
Nelson  (1995),  Hsu  (1996),  and  Kim  and  Nelson  (2006) 
present  methods  based  on  the  classical  statistical  model 
adopting  a  frequentist  view.  On  the  other  hand,  the  value  of 
information  (VIP)  approach,  exemplified  by  Berger  (1985), 
Bernardo  and  Smith  (1994),  Gupta  and  Miescke  (1996), 
Chick  and  Inoue  (2001),  and  Chen  et  al.  (2000)  uses  a 
Bayesian  framework  for  managing  the  trade-off  between 
the  consequences  of  an  immediate  decision  and  the  cost  of 
additional  sampling. 

The  overwhelming  majority  of  the  R&S  research  focuses 
on  a  single  unconstrained  performance  measure  (Kim  and 
Nelson  2006).  Most  applications,  however,  either  have  mul¬ 
tiple  performance  measures  (Butler,  Morrice,  and  Mullarkey 
2001)  or  face  some  constraints  on  the  primary  performance 
measure  (Andradottir,  Goldsman,  and  Kim  2005).  In  a 
manufacturing  setting,  for  example,  one  might  be  simulta¬ 
neously  interested  in  selecting  the  scheduling  policy  that 
maximizes  throughput,  with  a  physical  limit  on  the  work 
in  process  or  a  service  constraint  on  lead  times.  While 
some  of  these  measures  might  be  correlated,  others  might 
simply  conflict  with  one  another.  The  literature  dealing  with 
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multiple  performance  measures  or  stochastic  constraints  has 
been  sparse. 

In  the  case  with  multiple  performance  measures,  only 
a  partial  ordering  may  be  possible,  i.e.,  these  populations 
may  not  be  ordered  in  terms  of  a  vector-valued  parameter. 
In  such  cases,  one  typically  defines  a  real-valued  function  of 
the  parameters  and  use  this  function  to  rank  the  populations. 
Gupta  et  al.  (1973)  propose  several  such  scalar  functions. 
Andijani  (1998)  uses  the  analytic  hierarchical  process  (AHP) 
and  Butler  et  al.  (2001)  use  the  multiple  attribute  utility 
theory  (MAUT)  to  construct  a  scalar  function  by  assigning 
weights  to  the  performance  measures.  Kim  and  Lin  (2001) 
use  a  max-min  approach  whereby  they  maximize  the  poorest 
performing  criterion.  Lee  et  al.  (2006)  and  Lee  et  al.  (2007) 
seek  to  identify  a  Pareto  set  of  non-dominated  designs. 

R&S  in  the  presence  of  stochastic  constraints  has  re¬ 
cently  started  receiving  some  attention.  In  a  setting  where 
R&S  is  based  on  a  primary  performance  measure  subject 
to  the  feasibility  of  a  (possibly  correlated)  secondary  per¬ 
formance  measure,  Andradottir  et  al.  (2005)  propose  a 
two-phase  approach  :  phase  I  identifies  feasible  systems 
while  phase  II  selects  the  best  among  them.  With  the  ob¬ 
jective  of  accelerating  the  first  phase,  Batur  and  Kim  (2005) 
have  recently  introduced  procedures  for  finding  a  set  of  fea¬ 
sible  or  near-feasible  systems.  The  present  article  also  is 
concerned  with  efficient  feasibility  determination.  Unlike 
the  previous  work  that  is  based  on  the  indifference  zone 
perspective,  we  use  large  deviations  theory  to  minimize  the 
expected  number  of  incorrect  determinations.  We  charac¬ 
terize  the  optimal  employment  of  the  simulation  budget  to 
sampling  from  each  system,  and  present  a  stochastic  ap¬ 
proximation  algorithm  that  yields  a  budget  allocation  that 
provably  converges  to  the  optimal  one.  It  is  worth  pointing 
out  that  we  are  not  the  first  ones  using  the  large  deviations 
framework  for  R&S.  Glynn  and  Juneja  (2004)  have  used 
this  framework  not  only  to  determine  an  optimal  budget 
allocation  strategy,  but  to  also  demonstrate  the  deteriorating 
effect  of  violating  the  normality  assumptions,  which  are 
very  commonly  used  in  the  R&S  literature.  Our  notation 
closely  follows  theirs.  Like  them,  we  characterize  the  op¬ 
timal  allocation  of  the  computational  budget.  However,  we 
also  provide  an  algorithm  to  achieve  the  optimal  sampling 
allocations. 

The  remainder  of  the  paper  is  organized  as  follows. 
Section  2  introduces  the  necessary  notation  along  with  some 
preliminary  ideas  from  the  large  deviations  theory.  Section 
3  introduces  our  algorithm  for  feasibility  determination.  A 
numerical  illustration  is  presented  in  Section  4.  Section  5 
concludes  the  paper. 

2  OPTIMAL  ASYMPTOTIC  ALLOCATION 

We  consider  r  systems,  each  with  unknown  performance 
measure  , . . . ,  G  R.  Given  a  set  T  =  [7,°°),  we  wish 


to  employ  Monte  Carlo  simulation  to  determine  for  each 
system  i  whether  /r,  G  T  or  not.  We  assume  the  simulationist 
can  obtain  i.i.d.  replicates  of  the  random  variable  A,  to 
form  unbiased  estimators  A;(n,)  =  nJ^Ylk=\^i,k  of  jXi,  for 
i=  l,...,r. 

The  simulation  budget  n  is  allocated  in  order  to  mini¬ 
mize  the  expected  number  of  incorrect  determinations.  If, 
without  loss  of  generality,  we  assume  that  /Ti , . . . ,  /Tq  G  T  and 
llu+\ , . . .  ,IJ.r  then  the  feasibility  determination  problem 
is  to 

min  gn{pu---,Pr) 

Pi 

where 

g„(pi,...,Pr)  =  ^r)-|-  Y,  PiMPin)  &Y), 

and  ^  =  {pi  >  0 :  p,  <  1 }.  The  p;’s  represent  the  fraction 
of  the  simulationist’s  budget  that  is  allocated  to  sampling 
from  each  system,  where  for  simplicity  we  assume  that  each 
system  has  the  same  per-sample  cost. 

Our  contribution  is  two-fold;  (i)  we  characterize  frac¬ 
tional  allocations  p\,...,p*  that  are  optimal  (in  log  scale) 
as  n  ^  oo;  and  (ii),  we  provide  an  easily  implementable 
algorithm,  rooted  in  stochastic  approximation  theory,  that 
results  in  sampling  allocations  that  provably  achieve  the 
same  performance  (in  log  scale)  as  the  optimal  allocations 
in  the  limit  as  n  — >  oo. 

Let  /,(•)  be  the  large  deviations  rate  function 

/;(x)  =  sup{0x-logM,(0)},  (1) 

eeR 

where 

M,(0)  =£exp(0A;) 

is  the  moment  generating  function,  for  each  system  i  = 
Our  main  assumptions  are  that: 

AL  The  performance  measures  do  not  lie  exactly  at 
the  boundary;  /r,-  ^  7  for  all  i— 

A2.  For  all  systems,  Mi{d)  <°°  for  9  G  R. 

Assumption  Al  ensures  that  the  rate  functions  /,(•)  evalu¬ 
ated  at  the  boundary  point  7  are  positive,  so  that  no  system 
requires  all  the  simulation  budget.  Assumption  A2  is  stan¬ 
dard  in  situations  where  the  underlying  distributions  are 
light  tailed,  and  implies  that  for  each  system  there  exists 
9*  G  R  (unique)  such  that  7,(7)  =  9*Y  —  Mi{9*).  It  can 
be  shown  that  9*  <0  for  i  G  and  9*  >  0  for 

iG{a+l,...,r}. 

We  are  now  ready  to  state  the  main  results  of  this 
section. 
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Proposition  1  Suppose  Assumption  A2  holds.  Then  For  any  other  {p\,...,pr)  G  ^  we  have 


lim  -\oggn{p\,...,Pr)  = -minp,/;(7).  (2) 

n^oo  fi  I 

Proof:  The  proof  consists  of  two  parts.  For  Pi>0 
and  9  scalar 

-  \ogE  (exp(n0l,(p,n)))  =  p,  logM;(0 //?,). 
n 

Moreover, 


lim  inf  -\oggn{p\,...,Pr)  >  - - 

n^oon 

Proof:  The  first  statement  follows  immediately  by 

evaluating  Eq.  (2)  with  p*.  For  the  second  statement 
we  need  to  show  that  p\,...,p*  maximize  the  exponential 
decay  rate  of  gn(-),  given  in  Eq.(2).  Hence  we  consider  the 
problem  maXpg_^min,  p,7,(7),  which  is  equivalent  to  the 
linear  problem 


sup{70  -  p;logM,(0/p,)}  =  Pi  sup{y0/pi  -  logM;(0/p,)} 
eeJi  6eR 

=  pMr)- 

Hence,  the  G^tner-Ellis  Theorem  implies  that 
l/nlogP{Xi{npj)  ^  E)  ^  -Pikij)  for  i  = 
and  1  / n\ogP{Xi{npi)  G  E)  ^  —Pikij)  for  i  =  a  +  1 , . . . , r, 
as  «  — !■  oo. 

In  the  second  part  we  obtain  lower  and  upper 
bounds  for  (1/n)  logg„(pi , . . .  ,pr)-  Suppose  that  p,*/,*  = 
min,p,/;(7).  Then,  since  g„(pi, ... ,p^)  >  P(l,-(np;.)  ^ 
E)  if  i*  G  and  gn{p\,- ■  ■  ,Pr)  >  P{Xi*{npi*)  G 

E)  if  r  G  {a  +  1, . . .  ,r},  we  obtain  the  lower  bound 
liminf„^oc(l/n)logg„(pi,...,pr)  >  -Pi-h*(.7)-  For  the  up¬ 
per  bound, 

-\ogg„{pi,...,Pr)  <  -max{log(rP(li(npi)  ^E) ),..., 
n  n 

\og{rP{Xa{npa)  ^T)),\og{rP{Xi{npa+i)  GE)),..., 

\og{rP{Xi{npr)  G  E))}. 

Taking  limits, 

lim  sup  -\oggn{p\,...,Pr)  <  -pphif)- 


maxTj 

St  77  —  p,7,(7)  <0  for  i=  1 , . . . ,  r 

Ep;<1 

j=i 

p;  >  0  for  /  =  1 , . . . ,  r 
The  dual  of  this  problem  is 

min;r 

St  K  —  qdiiY)  >  0  for  i  =  1 , . . . ,  r 

E  7;  >  1 

j=i 

n,qi>Q  for  i—  l,...,r 

The  variables  pi  =  qi  =  ir' (y) (y),  r]  =  n  = 
I  {y)  are  primal  and  dual  feasible  and  r\  =  K, 

so  that  p*  is  primal  optimal  by  weak  duality,  (g) 

We  view  p*  as  the  optimal  allocation  scheme,  meaning 
that  no  other  allocation  achieves  a  higher  exponential  decay 
rate  for  g„{-)  as  the  sampling  budget  goes  to  infinity. 

Example  1  Suppose  for  every  system  we  have  Xi  ~ 
N{pi^af),  with  af  >  0.  For  each  system  i  (see  Exercise 
2.2.23  of  Dembo  and  Zeitouni  1998), 


g 

The  best  achievable  exponential  decay  rate  for  g„(-)  is 
described  next. 

Proposition  2  Suppose  assumptions  A1  and  A2 
hold.  Then 


lim  - 

n^oo  fi 


^Oggn{P*l,- 


,Pr) 


where 


* 


Pi  = 


if\y) 


(3) 


h{y) 


^  f Pi- y 


2 


Hence  Proposition  2,  specialized  to  the  Normal  case,  be¬ 
comes 


l!j=i  (yj/{Pj  -yV 


In  other  words,  when  the  underlying  distributions  are  nor¬ 
mal,  each  system  is  sampled  in  proportion  to  the  square 
of  the  distance  of  its  scaled  performance  measure  to  the 
boundary  point  7. 

Example  2  Assume  Xi  ^  Bernoulli{jii),  0  <  p,  <  1, 
for  i  =  1 , . . . ,  r,  and  let  0  <  y  <  1.  The  rate  functions  are 
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(Exercise  2.2.23  of  Dembo  and  Zeitouni  1998), 

so  that  Proposition  2  results  in 


2.  Generate  a  replicate  a  from  the  p.m.f. 
If;^/ZUlll^ori=\,...,r. 

3.  Update  sample  sizes:  Aa,n+i  =  ^a.n  +  1,  and 

“  ^i.n  for  i  f  cc. 

4.  Generate  a  replicate  from  system  a,  (say) 

5.  Update  Xa.„+i,  Mum+i,  Oa,n+u  and  Ia,„+i: 


7log 

i- 

)  -f  (l-7)log 

—  ) 
-Pj 

-1 

+-1 

7log( 

-1 

-f  (l-7)log 

+)] 

1  • 


(5) 


While  Proposition  2  characterizes  the  optimal  sampling 
allocation,  the  rate  functions  Ify)  are  generally  unknown 
and  need  to  be  estimated.  However,  estimating  each  rate 
function  involves  solving  a  root  finding  problem,  as  can  be 
seen  in  Eq.  (1)  by  solving  for  the  derivative  with  respect  to 
9  equal  to  zero.  This  suggests  that  the  computational  cost 
of  sequentially  estimating  the  rate  functions  dominates  all 
other  costs  when  running  the  simulation,  and  therefore  is 
not  a  promising  course  of  action.  In  the  next  section  we 
present  a  stochastic  approximation  algorithm  that  overcomes 
these  issues,  and  leads  to  fractional  allocations  that  converge 
almost  surely  to  the  optimal  p*  allocations. 


3  STOCHASTIC  APPROXIMATION  ALGORITHM 


Initially  we  warm-up  each  system  with  n^  >  0  replicates 
and  compute  the  sample  averages  Xifi  =  Mq  *  Let 

f  /  1 

/,, 0  =  sup  <  07- log  —T  exp{eXi^k) 

=  Oi fiY- log  (Mifi) , 


where  0,  o  satishes  the  root  problem 

_  lHUXi^k^y^piOifiXLk) 

i;iiexp(0,pA,-,,)  ’ 

and  the  sample  moment  generating  functions  M,  o  evaluated 
at  0j',o  are  given  by 


2  no 

^i.o  =  —  y"  exp(0,pA,  ,(-). 

no  it) 

In  order  for  the  algorithm  to  get  properly  initialized,  if 
/,p  <  0  we  set  7,  0  =  e,  for  e  a  small  positive  constant. 

Finally,  we  initialize  the  sample  sizes,  A, p  =  no  for 
i  =  l,...,r. 

Algorithm. 


i 


Mr 


,71+1  — ^o:,«  “1“  ^  i.^^P{^(X,n^a,Xa^  ^a,n)i 

'^a.Ti+l 


0ce,n+l  —  0a,n 

-  -  7) exp  {0a,„{Xa^l^  -  7))  , 


and 


^a,7i+l  —  ^a,7i+l  7-log(Ma.„+i), 

if  0a,«+i7-log(Ma,„+i)  >  0,  and/a, „+i  =  miny/p„ 
otherwise.  (This  ensures  that  Ia,n+i  >  0  always.) 
For  i  cc,  set  —Xi  ^i,  Mj  fj^i  =Mi  fj,  0/,«+i  = 
^i,n  7  and  Ii^n+\  —  f,n  ■ 

6.  Increase  n^n+l  and  go  back  to  2. 

To  ensure  that  each  system  is  sampled  inhnitely  often, 
let  (v„)  be  an  increasing  sequence  such  that  v„  °°  and 
n^'lLi7(Vi  <  n)  — >  0,  where  Jf)  is  the  indicator  function. 
We  sample  from  all  r  systems  at  stages  Vi ,  V2, . . .,  and  update 
the  parameters  according  to  steps  3  and  5  of  the  algorithm. 

This  algorithm  provides  fractional  allocations  that  con¬ 
verge  a.s.  to  the  optimal  allocations,  as  the  next  theorem 
demonstrates. 

Theorem  1  Under  assumptions  A1  and  A2, 


almost  surely  as  n  — >  00. 

Proof  outline:  Suppose  momentarily  that  9i  „  —f  9* 
a.s.  Let  pi  „  =  Xiftjn,  so  that  step  3  of  the  algorithm  can 
be  expressed  as  p,>+i  =  -f  (7(a„  =  /)  -  pi,n)/{n+  1), 
where  a„  is  the  nth  replicate  of  a  generated  in  step  2  of  the 
algorithm,  and  7(  )  is  the  indicator  function.  The  recursion 
for  pi^n+\  can  be  re-written  as 


Pi,n-\-\  —  Pi,n  T  _  ,  ,  (Pj  Pi,n)  “f  _  ,  ,  —  0  9i,n) 


n-f  1 


n  -f  1 

+  — pr+n-Pi*):  (6) 
n-\-  i 


1.  Initialize  n  =  0. 
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where  qi^„  =  lij; /YJj=\^j)i-  Since  = 

=  0,...,n]  =  the 

sequence  {7(a„  =  i)  —  qi^„}  is  a  martingale  difference 
with  respect  to  the  sequence  of  a-algebras  gener¬ 
ated  by  =  0,. 

Additionally, 


The  next  result  shows  that  the  expected  number  of 
incorrectly  determined  systems  produced  by  the  stochastic 
approximation  algorithm  approaches  0  at  the  best  possible 
rate. 

Theorem  2  Suppose  that  assumptions  A1  and  A2 
hold.  Then, 


l^i.k  Pi  I  ^  Pi  I  >1/(1+^)) 

r  ^+1  ~k  ^+1 

almost  surely.  The  set  {T.kJi\^i k  —  pl\  >  l/(l+^))/(^  + 

1)  =  -}  c  >  1/(1  +^))/(^+ 1)  =  -}, 

and  P{i:kJ{\qi,k-p*\  >  1/(1  +  f))/(^  +  1)  = -)  =  0  for 
each  £  since  Qi  n  9*  a.s.  implies  qi^„  p*  a.s.  Therefore 
PiLkJ{\qi,k-  P*\  >  1/(1  +^))/(^+  1)  =  °°)  =  0,  and  we 
conclude  from  Eq.  (7),  that  Y.k  \qi,k  —  P*i\/{k+  1)  <  oo  a.s. 

Hence,  all  the  assumptions  of  Theorem  5.2.1  ofKushner 
and  Yin  (2003)  are  satisfied.  Moreover,  the  ordinary  dif¬ 
ferential  equations  =  p*  —  pt,  i  =  I, . . .  ,r,  have  a  unique 
globally  asymptotically  stable  point,  p*,  so  that  Theorem 
5.2.1  of  Kushner  and  Yin  (2003)  results  in  „  ^  p*  a.s. 
(The  special  instances  of  n  =  Vk,k=  1,2,...,  can  be  handled 
by  introducing  an  extra  term  in  Eq.  (6),  and  do  not  affect 
the  asymptotic  behavior  of  p,;„.) 

Regarding  the  asymptotic  behavior  of  letg,(0,'  „)  = 
^[(Y,' —  7)exp(0,'.„(Y,' —  7))]  and  re-write  the  recursion  for 
9i,n  as 


9a.n+l  —  9a, n  +  ^ -  X  8a{9a,n) 

-  iXa,la  -  7)exp  {9a,n{Xa,l,  -  7)) )  - 

/  ^a,n+l 

The  sequence  {8a{9a,n)  -  {Xa,la  “ 

7)  exp  (0c(,n(Ya  —  7))}  is  a  martingale  difference 

with  respect  to  the  sequence  of  a-algebras  generated 
by  =  0,...,n}.  Since 

Xin  °°  a.s..  Theorem  5.2.1  of  Kushner  and  Yin  (2003) 
applies,  meaning  that  0,  „  converges  a.s.  to  the  unique 
globally  asymptotically  stable  point  of  the  ODE 

9't  =  -m), 

which  is  given  by  0/ ,  for  i=  1 , . . . ,  r.  (g) 

Remark.  P*iU{Y)  =  pYiiy)  suggests  that  step  2  of 

the  algorithm  could  be  replaced  by  a  =  arg min, {A,- 
An  argument  similar  to  the  one  employed  in  the  proof 
of  Theorem  1  shows  that  this  approach  also  results  in 
asymptotically  optimal  fractional  allocations. 


1 


log  ^p(Y,-„^r)+  £  p(Y,-„Gr) 

\!=1  i=a+l  J 


as  n  °o. 

Proof  outline:  The  proof  resembles  the  proof  of  Propo¬ 
sition  1.  Eirst  of  all,  the  random  variables  Y,  1, . . .  ,Y, 
are  i.i.d.  conditional  on  A,  „,  so  that 


1 

n 


log(£'exp(n0Y,-_„)) 


-logE 

n 


Mi 


where  the  expectation  is  with  respect  to  the  probability 
measure  of  the  random  variable  A,_„.  Assumptions  A1  and 
A2,  Theorem  1,  and  the  continuity  of  M,(-)  lead  to 


Mi  0 


Xi,n 


hn/n 


>Mi 


Q\P. 

Pi 


(8) 


a.s.  By  Jensen’s  inequality 


Mi\  0 


^i,n 


<E 


Mi  0 


^i,n 


SO  that 


p*  logM,-  (  )  <  lim  inf  -  log  {E  exp(n0Y,-  „)) .  (9) 

\p1J 


By  Eq.  (8),  given  e  >  0, 


Xin 


X,nl^ 


Mi[9—\  <Mi[- 


0  \ 


Pi 


a.s.  for  all  n  sufficiently  large,  which  results  in 

lim  sup  -  log  {E  exp(n0Y,;„))  <  {p*  +  e)  logM,-  (  -^  )  ■ 
n^oo  n  \Pi  J 


Sending  e  ^  0,  together  with  Eq.  (9),  yields 


-  log  {E exp(n0Y,;„))  ^  p*  logM; 


By  the  G^tner-Ellis  theorem,  we  conclude  that 
l/nlogP(Y;_„  ^  r)  ^  —Ppiiv)  for  i  —  l,...,a,  and 
l/nlogP(Y;_„  G  r)  ^  —Ppiiy)  for/  =  a-f  1, . . .  ,r,  as  n  ^  00. 
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The  rest  of  the  proof  is  similar  to  the  second  part  of  the 
proof  of  Proposition  1,  and  is  thus  omitted.  0 

4  NUMERICAL  EXAMPLES 

We  now  illustrate  our  stochastic  approximation  algorithm 
through  the  two  examples  introduced  in  Section  2,  where 
the  performance  metric  of  interest  is  distributed  according  to 
Normal  and  Bernoulli  distributions,  respectively.  In  addition 
to  providing  us  with  an  opportunity  to  examine  the  hnite- 
time  performance  of  our  algorithm,  these  two  examples  also 
illustrate  the  impact  of  deviations  from  the  assumption  of 
normality,  which  is  quite  common  in  the  R&S  literature. 
All  the  results  are  the  average  of  100  independent  mega¬ 
replications.  The  results,  which  reflect  the  allocation  of  the 
computing  budget  as  a  function  of  the  number  of  iterations 
of  the  algorithm,  lead  to  two  key  observations:  the  efficacy 
of  sequential  algorithms  and  the  importance  of  adequate 
initialization  of  such  algorithms.  The  efficacy  of  these 
algorithms  is  particularly  valuable  in  critical  settings,  i.e., 
those  settings  where  the  system  performance  is  close  to  the 
boundary  of  the  feasible  region. 

Example  3  Suppose  we  have  five  systems  with 
Xi^N{pi,af),  with  ^,  =  [9.51,9.45,9.40,9.55,9.60]  and 
of  =  1.  Let  the  boundary  point  J—  9.50.  As  a  starting 
point,  we  will  use  Eq.  (4),  and  sample  from  each  system 
in  proportion  to  the  square  of  the  distance  of  its  scaled 
performance  measure  to  y.  In  this  case,  the  optimal  budget 
allocation  for  the  five  systems  is  given  by  p*=  [0.9091, 
0.0364,  0.0091,  0.0364,  0.0091].  Figure  1  shows  that  con¬ 
vergence  is  rather  rapid,  that  is,  within  100000  iterations, 
the  optimal  allocation  of  the  computing  budget  is  achieved 
along  with  the  correct  classification  of  all  the  systems.  Note 
that  the  algorithm  has  been  initialized  with  only  100  samples 
from  each  system. 


Figure  1:  Convergence  of  the  budget  allocation  as  a  function 
of  log(iterations)  for  the  Normal  case  with  the  exact  rate 
function 


Given  the  fact  that  the  rate  functions  are  generally 
unknown  and  need  to  be  estimated,  we  next  deploy  our 
stochastic  approximation  algorithm  to  classify  five  systems 
withXi^N{pLi,af),  with  Pi=  [9.20,8.50,9.00,9.80,10.01] 
and  <jf  =  1.  Let  the  boundary  point  7=9.50.  In  this 
case,  the  optimal  budget  allocation  for  the  five  systems  is 
given  by  p*=  [0.3577,  0.0322,  0.1288,  0.3577,  0.1238]. 
Figure  2  shows  that  convergence  has  slowed  down  due  to 
the  effort  required  in  estimating  the  rate  functions;  that 
is,  it  now  takes  over  1000000  iterations  for  the  stochastic 
approximation  algorithm  to  achieve  the  optimal  allocation  of 
the  computing  budget  along  with  the  correct  classification  of 
all  the  systems.  Note  that  the  algorithm  has  been  initialized 
with  50000  samples  from  each  system. 


“■-•System  1 
-■-System  2 
-5  Systems 
System  4 
Systems 


Figure  2:  Convergence  of  the  budget  allocation  as  a  function 
of  log(iterations)  for  the  Normal  case  with  the  estimated 
rate  function 

Example  4  Assume  A,-  ~  Bernoulli{p.i),  with  Pi  = 
[0.92,0.85,0.90,0.98,0.88].  Let  the  boundary  point  y  — 
0.95.  Based  on  Eq.  (5),  the  optimal  budget  allocation  for  the 
five  systems  is  given  by  p*  =  [0.4492,  0.0618,  0.1878,  0.1927, 
0.1084].  Figure  3  shows  that,  within  10000  iterations,  our 
stochastic  approximation  algorithm  achieves  the  optimal 
allocation  of  the  computing  budget  along  with  the  correct 
classification  of  all  the  systems.  One  must  recall,  however, 
that  appropriate  initialization  of  the  algorithm  is  crucial. 
The  above  performance  is  obtained  after  collecting  50000 
observations  to  initialize  the  algorithm. 

5  CONCLUDING  REMARKS 

We  address  the  problem  of  feasibility  determination  in  a 
stochastic  setting.  More  specifically,  we  wish  to  determine 
whether  a  system  belongs  to  a  given  set  F  based  on  a  perfor¬ 
mance  measure  estimated  through  Monte  Carlo  simulation. 
Our  contribution  is  two-fold:  (i)  we  characterize  fractional 
allocations  that  are  asymptotically  optimal;  and  (ii)  we  pro¬ 
vide  an  easily  implementable  algorithm,  rooted  in  stochastic 
approximation  theory,  that  results  in  sampling  allocations 
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Figure  3;  Convergence  of  the  budget  allocation  as  a  function 
of  log(iterations)  for  the  Bernoulli  case  with  the  estimated 
rate  function 


that  provably  achieve  in  the  limit  the  same  performance  as 
the  optimal  allocations.  The  finite-time  behavior  of  the  algo¬ 
rithm,  illustrated  on  two  small  examples,  is  quite  promising. 
Current  research  is  aimed  at  generalizing  the  approach  to 
settings  where  the  feasible  region  is  cZ-dimensional  with 
d>\. 
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